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SUMMARY 


The quadratic isoparametric elements which embody the inverse square root 
singularity are used for calculating the stress intensity factors at tips of 
cracks. The strain singularity at a point or an edge is obtained in a simple 
manner by placing the mid-side nodes at quarter points in the vicinity of the 
crack tip or an edge. These elements are implemented in NASTRAN as dummy ele- 
ments. The method eliminates the use of special crack tip elements and in 
addition, these elements satisfy the constant strain and rigid body modes re- 
quired for convergence. 


INTRODUCTION 


In "Crack Tip Finite Elements are Unnecessary", Henshell and Shaw (ref. l) 
reported that the inverse of the Jacobian associated with the coordinate 
transformation becomes singular at a point when the mid-side nodes for two- 
dimensional eight-point quadrilateral elements are placed at quarter points. 
Interestingly enough, the same singularity was discovered independently by 
Barsoiim (ref. 2) for two-dimensional, as well as, three-dimensional quadratic 
isoparametric elements. It was then natural to investigate the order of the 
singularity and it was fotind that the singularity was precisely of the order 
one-half for the strains , a phenomenon encountered in linear fracture mechan- 
ics. This remarkable phenomenon completely eliminates the necessity of incor- 
porating special crack-tip elements (ref. 3> ^ and 5) and has additional ad- 
vantages over the special crack-tip elements, namely; it satisfies constant 
strain and rigid body modes. The special crack-tip elements were introduced 
in the literature to avoid the extremely fine grid mesh required in the vicin- 
ity of the crack and the cumbersome extrapolation needed when using regular 
finite elements (ref. 6 and 7). 

Advanced versions of NASTRAN (ref. 8), as well as some general purpose 
programs have such isoparametric elements. Hence, by Judicious choice of 
nodes, accurate crack-tip elements can be formulated and stress intensity 
factors for cracks and flaws can be computed. 

In this paper, after a brief review of the two and three-dimensional 
formulation, we discuss the implementation of the two-dimensional quadri- 
lateral. and three-dimensional brick elements as NASTRAN duimny user elements. 


I 



Lastly, test problems are done to assess the accuracy. Stress intensity 
factors are computed for a C-shaped specimen. The C-shaped fracture toughness 
method has been accepted by ASTM as a standard test for thick-walled cylinders. 


SYMBOLS 


(x,y) 

Cartesian coordinates 

(C,n) 

Curvilinear coordinates 

Xf ,y^ ,Ci 

Grid point coordinates 

Ni 

Shape function at grid point i 

u,v 

Cartesian displacements 

{e} 

Strain vector 

[J] 

Jacobian matrix 

{a} 

Stress vector 

[D] 

Stress-strain matrix 

[K] 

Element stiffness matrix 

E,G,v 

Elastic constants 

Kl.Kii 

Stress intensity factors 

r.e 

Local cylindrical coordinates 

{F}e 

Equivalent nodal forces 


THE TWO-DIMENSIONAL CASE 


Following the notation of Zienkiewicz (ref. 9) the eight node element in 
Cartesian coordinates (x,y) is formulated by mapping its geometry into the 
c\irvilinear space (C,n) of the normalized square (-1 1 C 1 1, -1 1 n £ l) by 
quadratic shape ftmctions of the ^Serendipity ’ family (ref. 9): 


X — Z N (5,r))xj^ , 

i=a 
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y = I Ni(c,n)yi » 

(1 

Ni = [(i+C5i)(l+firii) - (i-C^)(l+nni) - (l-n^)(i+CCi)] cf 

+ (l-5^)(l+nTii)(i-C?)Ti?/2 + (i-Ti2)(i+5Ci)(i-Ti|) 5?/2 , 

viiGrc IS *tiiG slisips function nt no^6 i whoso C&i*t©si&n snd cnrviXino&r co™ 
ordinates are (xi,yi) and (Ci,ni) respectively. The details of the shape 
fxinctions and the numbering sequence are given in figiire 1. The same shape 
fTonctions are used to interpolate the displacements within the element, hence 
the name isoparametric : 


u = 2 Ni(c,n)ui , 

i=l 


u 

V = E Ni(5,ti)vi . 
i=l 

The stiffness matrix is found in the usual way as follows: 


Substituting from equation (2) into equation (3) we have! 


1 • 


1 

1 • 

1 


1 

) 

> — [ • • • • • • ] ^ 

^i I 

J V* 

^i 1 

t 1 

( 1 


1 • 

1 • 

\ • i 


• 1 


where 


[B^] = 


3Ni 

0 

w 

0 

3Ni 

9y 

3Ni 

3Nj 

3x 

9y 
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By the rules of partial differentiation we obtain 


3Ni 


3Ni j 

3x 1 


\h j 


' = [J]-1 

3Nj 1 

j 

|3Ni j 

3y ^ 

1 

) 


( 6 ) 


where [j], the Jacobian matrix, by virtue of equation (l) is given by 

I I 


[J] = 


3x 



9C 

H 




sc 

3x 

3jr 


9n 

3n 



3Ni 

W 

3n 


I I 

^i yi 


I 


( 7 ) 


The stress components are given by 
cr-r 


{a} = 


= [D] {e} 


( 8 ) 


^xy 


where [D] is the stress-strain matrix and for the case of plane stress is 
given by 


[D] = 


E 


l-v2 


1 V 0 
V 1 0 

0 0 (l-v)/2 


( 9 ) 


The element stiffness matrix is then: 


[K] = /:/: [B]^[D][b] det|j|dCdn . 


( 10 ) 


The integration in equation (lO) is done numerically by nine-point Gaussian 
quadrature as explained in reference 9- 


THE CRACK TIP ELEMENT 


It is clear from equations (U) and (6) that we rieed the inverse of 
Jacobian matrix [j] before the strains can be computed. Hence, whenever the 
inverse of [j] is singular or, equivalently, the determinant of [j] is zero, 
the strains and stresses become singular. This is simply accomplished by 
placing the mid-side nodes (e.g., nodes 5 and 8 of figure l) at quarter points 
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from node 1 in Cartesian coordinates. 


This can be illustrated by investigating the singularity along line 1-2 
(n=-l) of fig\ire 1. Evaluating the shape functions given in figure 1 at n=-l. 


we have the transformation 

X = -1/2 5(1-5) + 1/2 5(1+5) X2 + (1-5^) X5 . (11) 

Choosing xq^=0, X 2=L and the quarter point x^=L/k, equation (ll) becomes 

X = 1/2 5(1+5)L + (1-5^) L/h (12) 

Solving for 5 we have 

5 = -1 + 2/5^ . (13) 

In this case the reduced Jacobian becomes 

|| = i (IH) = ^ . (IM 


Equation (lJ+) clearly indicates the singularity for the inverse of the Jacobian 
at x=0, 5=-l* The order of singularity can be obtained from the displacement 
along line 1-2 (fig. l). From equation (2) we have 

u(5,-l) = -1/2 5(l-5)ui + 1/2 5(1+5 )u 2 + (l-5^)u5 , 
and writing in terms of x from equation (13) we have 

u = -l/2(-l+2*^)(2-2»^)ui + l/2(-l+2>^)(2>^)u2+lt(*^ - §)u^ . (l5) 

Differentiating equation (15) we obtain the strains in the x-direction: 

tx = ^ ri“2*t4 - ^1U5 , (16) 

indicating the singularity of order one-half (^ )> precisely the singularity 

needed for crack problems. It can be seen that equation (16) also incorporates 
constant strain terms. 

We have only investigated the singularity at node 1 along line 1-2 of fig- 
ure 1. However, the singularity at node 1 along any other ray emanating from 
node 1 is weaker than one-half. The singularity of order one-half can simply 
be achieved by collapsing grid points 1, 1+ and 8 and placing grid points 5 and 
T at the quarter points in Cartesian coordinates as shown in figure 2. With- 
out loss of generality we take the Cartesian coordinates as shown. Using equa- 
tions (1) and (7) it can be shown that 

det. |j| = 1/16 (l+5)^sina , (iT ) 

which vanishes for 5=-l for all n (i.e., along any ray from node l). The 




di splEc 611161113 , using polar coordinates (x = rcosG, y = rsin6) is given liy; 


U = 2rl/2 

cosl/2(a/2) 




-i/i+(i+n)u 2 + l/ 2 (i-r))u 5 + i/ 2 (i+n)uYl -i/ltn(i-n )u 2 +i/i^n(i+n)u- 


+1/2(1 


with 


-n^)u6 


_ tan(e- g/2) 
tan (a/2) 


(18) 


(19) 


E(juation (18) indicates that the strains will have the necessary singularity 
of order one-half (^/i) at the crack tip. 

The stress intensity factors Kj and Kjj are computed at the quarter points 
using the Westergaard near field displacements which, for plane stress, are 
given by (ref. 10 ). 


Vt r. 1/2 

U =1 (^) COS0/2 

G 27 T 


+ cos^e/2 


I 


+ sin20/2l + ^^^(~)l/2sin0/2 - 
[l+v J G ] 


l+v 


V 

l+v 


= (— )^'^^sin0/2 f_^ - cos2e/2l + ^^l(_2L)^^^cos0/2’ T- Iz. 

TT 2 ir [^ 1 +v . J G 2 ir L 1 + 

] 


( 20 ) 


+ sin^0/2 
Solving for Kj, Kjj we have: 


u COS0/2 [- + cos^0/2] + V sin0/2 + cos^0/2] 


Kj = G(^)1/2 
-L 


l+v 


l+v 


- cos2e/2][y|;7 - cos2e/2] 


2ir 1/2 u sin0/2 - v cos0/2 
KlI = G(--) 


( 21 ) 


l+v 


THE THREE-DIMENSIONAL CASE 

The three-dimensional twenty-point isoparametric quadratic 'Brick' element 
is formulated in much the same way, by mapping the geometry into curvilinear 
space (C,n,;) of a normalized cube (-1 1 5,n,C 1 l) by the quadratic shape 
function (ref. 9), 
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X 


9 


20 

= E Ni(C,ri»C)xi 
i=l 

y = E° Ni(c,n,c)yi , 
i=l 

z — Z Nj^ ( ^ }T1 } C ) t 
i=l 

% = l/8(i+C5i)(l+nni)(l+?Ci)(CCi + + CCi -2) cf (22) 

+ lA(i-52)(l+ntii)(i+C?i)(l-c2) 

+ iA(i-n2)(i+«i)(i+cCi)(i-n2) 

+ lA(l-c2)(l+|5j)(l+,1Tli)(l-c2) , 

where is the shape fvinction at node i (i=l to 20) whose Cartesian and cvirvi- 
linear coordinates are (xj, yi, zi) and (Ci, rii» ?i) respectively. It should be 
noted that the shape function given in equation (22) is obtained by superposi- 
tion of those given in reference 9* The geometry of the unit cube and the ntim- 
bering sequence, as suggested by reference 11, is shown in figure 3. 

For the isoparametric fomulation and displacements are given by 
20 

u = Ni(5,b,?)ui , 

20 

■V = %(C,n,C)vi , (23) 

20 

V “ Ni(C,n,c)wi . 

The rest of the analysis follows in a similar fashion that given for the two- 
dimensional case with appropriate augmentation to the three-dimensional 
quantities. 

The singularity element is obtained by collapsing one face, 2376, and 
placing the midside nodes 9» 13, 11, 15 of figure 3 at quarter points. The 
singular element is shown in figvure 4, in Cartesian coordinates. Since the 
elements are isoparametric they automatically satisfy inter-element compati- 
bility and continuity in their regular or singular forms. It should be noted 
that the displacements are not singular. Further, it is easily shown that 
I Ni=l. Hence, by theorems given in reference 9 the elements satisfy the con- 
stant strain and rigid body modes. The above conditions are necessary for the 
'patch test' mentioned in reference 2. 
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NASTRAU IMPLEMENTATION 


The isoparametric quadratic quadrilateral and brick elements have been im- 
plemented using the NASTRAN dummy user element facility as outlined in Section 
6 . 8.5 of reference 12. This involved coding element stiffness and stress data 
recovery subroutines using the analysis outlined above and relinking the 
affected NASTRAN Links. Additional modifications were required to some of the 
Output File Processor COFP) routines. These changes are detailed below. 

The quadrilateral element was implemented as a DUMl element. Figure 5 
shows the formats for the ADUMl, CDUMl and PDUMl cards. KDUMl, the element 
stiffness matrix subroutine, obtains material and grid point information from 
the element connection and property table (ECPT) and builds the matrices re- 
quired to perform the integration in equation (lO). The integration is per- 
formed numerically using compound 3-point Gaussian quadrature as explained in 
reference 9* This results in 9 evaluations of the integrand. Once the l6byl6 
stiffness matrix is complete, the appropriate 2 by 2 submatrices corresponding 
to the given pivot point are entered into the upper left of the 6 by 6 sub- 
matrices required by SMAIB, the stiffness matrix insertion subroutine. SMAIB 
is called 8 times for each pivot point. The time for element stiffness genera- 
tion is lU seconds per element on an IBM 36 O model ItU. 

NASTRAN stress data recovery is accomplished in two phases. During phase 
I, SDUMll calculates [D][b] from equations ( 5 ) and (9) for each grid point and 
passes the resultant 2k by 16 matrix to SDUM12 for final stress calculations. 
SDUMll also checks for singularities in the inverse of the Jacobian (eq.(7)) 
and flags those grid points which have a singxxlarity . Information passed to 
SDUM12 from SDUMll includes the element id, grid point numbers, grid point 
singularity flag, coordinates of the eight grid points and the material con- 
stants E, G and v. SDUM12, phase II of the stress recovery, locates the dis- 
placements associated with a given element and multiplies [D][B] times these 
displacements (eq. (It) and (8)) to give the stress components at each of the 
eight grid points. Grid point flags are checked for singular grid points and, 
if singularities exist. Mode I and Mode II stress intensity factors are calcu- 
lated using equation (21). These stress intensity factors at the quarter 
points are output at the corner nodes of the collapsed side while the corres- 
ponding mid-side node stress output is set to zero. The point ids of the 
singular corner nodes are negated and the mid— side id is set to overflow the 
integer field specification, thus flagging the point with asterisks. 

OFF has been modified to output the eight sets of stress components for 
each element. These modifications were implemented by adding heading formats 
to OFPIA and changing the appropriate pointers and format specifications in 
OFPIBD, 0FP5BD and OFIPBD. Although the ADUM cards allow sufficient flexi- 
bility to implement the element stiffness subroutine, changes were required to 
GPTAl which describes the connection and property characteristics of each ele- 
ment (see Section 2. 5.2.1 in ref. 12). These changes were required to handle 
the expanded stress requirements. The number of words SDR2 passes from phase 
I to phase II was changed from 100 to 1+30 while the count of words SDR2 out- 
puts for real stresses was increased from 10 to 33. 
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The implementation of the twenty node brick element was not as straight- 
forward as that of the quadrilateral element. The brick was implemented as a 
DUM2 element. Figiare 6 shows formats for the ADUM2, CDUM2 and PDUM2 cards. 
DIMENSION changes were required in TAIA and TAIB since NASTRAN assumes a maxi- 
mum of 10 grid points per element. KDUM2 was initially implemented similarly 
to KDUMl. However, in this case, the integration using 3-point Gaussian quad- 
rature requires 2? integrand evaluations with a stiffness matrix of order 60. 
The size of the KDUM2 subroutine necessitated a change to the overlay structure 
of LINKS, placing KDUM2 in its own overlay segment. Also, since NASTRAN calls 
the stiffness routines based upon the pivot point concept, the same brick ele- 
ment stiffness matrix is built twenty times in an analysis. This technique 
results in a 20 minute stay in SMAl for a one element problem on an IBM 360 
model 44. Changes are being made to KDUM2 to build each element stiffness 
matrix once and save it on auxiliary storage. When a request for an element 
stiffness matrix is made, KDUM2 will check auxiliary storage for a copy of the 
matrix. If it is not there, ICDUM2 will build the matrix and add it to the 
file. If it is there, the matrix will be retrieved and not recalculated. Us- 
ing this procedure, stiffness matrix generation should not exceed 2 minutes 
per element . 


Stress data recovery is also non-standard for the brick element. Due to 
the size of the arrays used in stress recovery for the brick, phase I has a 
limited function of assembling arrays dependent on parameters not available to 
phase II. The majority of calculations for stress recovery are accomplished 
during phase II, saving storage but increasing stress recovery times. Quanti- 
ties passed from phase I to phase II are the stress-strain matrix [D], the 
element id, grid point ids, grid point coordinates and the material constants 
E, G and v. Phase II locates the displacements, calculates the stress compo- 
nents for each grid point, checks for singular Jacobians and calculates stress 
intensity factors as required. Stress intensity factors are displayed in a 
manner similar to that employed for the quadrilateral element. 

OFP has been modified to output the twenty sets of stress components for 
each element. A DIMENSION change was also required in OFP to allow twenty 
grid points per element. OFPIA, OFPIBD, 0FP5BD and OFIPBD were updated to 
produce the required heading and output formats. GPTAl was changed to increase 
the nxjmber of words to 120 that SDR2 passes from phase I to phase II. The 
count of words SDR2 outputs for real stresses was increased to l4l. 

Both element implementations were checked independent of NASTRAN via dummy 
driver routines for SMAl and SDR2. The coding for the element stiffness sub- 
routines was verified by multiplying the stiffness matrix for one element times 
the known displacements for a uniform stress field. The resultant nodal load- 
ing was compared to that found analytically from (for the two-dimensional case) 



{a} d? dn 


(24) 


Figure 7 illustrates equivalent nodal forces for Oy=l on the normalized square 
(~1 i x,y il). The stress recovery coding was checked by passing the known 
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displacements for a iiniform stress field to the stress recovery subroutines 
and observing the constant stress results. 

Implementation of the quadrilateral element was complete within a week of 
the completion of the analysis while the brick element implementation took 2 
weeks because of the length of coding, overlay changes and NASTRAN changes to 
support more than 10 grid points. 


NUMERICAL RESULTS 


To assess the accuracy of the method, three test problems with fairly a 
coarse grid (68 elements and approximately 239 grid points) were run. The 
problems are the single edge crack, double edge crack and center crack. These 
three problems can be done in a single run as subcases with different single 
point constraints as shown in figures 8a, 8b, 8c. The solutions of above 
problems, by various methods, have been well documented (ref. 12). It was 
foxmd that the finite element solutions were accurate to within 2 - 3 %. Graphic- 
s^ly » this is illustrated in figure 9 for the double edge crack by using the 
Westergaard near field solution (ref. 10) for Oy. 

ASTM has stringent requirements for the size of specimens for fractxire 
toughness testing. However, in many applications of thick-walled cylinders 
these requirements are not easily met. The C-shaped specimen, which is easily 
obtained from thick-walled cylinders, was suggested (ref. 13) and is now 
accepted as a standard test for such cylindrical material. The stress intens- 
ity factors for such a section, shown in figure 10, were computed for different 
crack lengths and the finite element results, experimental results, and the 
collocation results of reference 13 are shown in figure 11. It is seen that 
remarkable agreement is obtained with just U8 elements and 171 grid points. 

The results have also been compared with those in reference l6 and similar 
correspondence was observed. 


CONCLUSIONS 


Quadratic isoparametric elements have been used to form a singular ele- 
ment for fracture mechanics analysis. These elements provide excellent results 
even with coarse grids as long as the singular elements strictly conform to the 
geometries of figures' 2 and 1+. The elements have been successfully implemented 
on NASTRAN Level 15*0 as dummy user elements. 
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i/U(i+c)(i+n)(?+n-i) 



and Nvunbering Sequence for 







Figure 5. Bulk Data Cards for Quadrilateral Element. 



Figure 6. Bulk Data Cards for Brick Element. 
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Figure 8. - Test Problem: Shaded Area Analysed by Finite Elements. 





2.5 



e- DEGREES 


Figure 9« - Comparison of NASTRAM Results With the Theoretical 
Solution for Double-Edge Crack. 
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